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Abstract 

It  is  known  that  small  relative  perturbations  in  the  entries  of  a  bidiagonal  matrix 
only  cause  small  relative  perturbations  in  its  singular  values,  independent  of  the  values 
of  the  matrix  entries.  In  this  paper  we  show  that  a  matrix  has  this  property  if  and  only 
if  its  associated  bipartite  graph  is  acyclic.  We  also  show  how  to  compute  the  singular 
values  of  such  a  matrix  to  high  relative  accuracy.  The  same  algorithm  can  compute 
eigenvalues  of  symmetric  acyclic  matrices  with  tiny  componentwise  relative  backward 
error.  This  class  includes  tridiagonal  matrices,  arrow  matrices,  and  exponentially  many 
others. 


1  Introduction 

\ 

In  [9]  it  was  shown  that  small  relative  pertu£bations  in  the  entries  of  a  bidiagonal  matrix 
B  only  cause  small  relative  perturbations  in  its  singular  values.  This  is  true  independent 
of  the  values  of  the  nonzero  entries  of  B.  This  property  justifies  trying  to  compute  the 
singular  values  of  B  to  high  relative  accuracy,  and  is  essential  to  the  error  analyses  of  the 
corresponding  algorithms  [9]. 

Since  this  attractive  property  of  bidiagonal  matrices  is  independent  of  the  values  of  the 
nonzero  entries,  it  is  really  just  a  function  of  the  sparsity  pattern  of  bidiagonal  matrices. 
In  this  paper  we  completely  characterize  those  sparsity  patterns  with  the  property  that 
independent  of  the  values  of  the  nonzero  entries,  small  relative  perturbations  of  the  matrix 
entries  only  cause  small  relative  perturbations  of  the  singular  values.  The  characterization 

‘The  author  was  supported  by  NSF  grant  ASC- 9005933  and  DARPA  grant  DAAL03-91-C-0047  via  a 
>ubconiiacl  from  the  University  of  Tennessee.  This  work  was  performed  during  a  visit  to  the  Institute  for 
.Malhenialic.s  and  its  Applications  at  the  University  of  Minnesota. 

'  riiv  author  also  acknowledges  the  Institute  for  Mathematics  and  its  Applications  at  the  University  of 
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is  simple:  a  sparsity  pattern  has  this  property  if  and  only  if  its  associated  bipartite  graph 
is  acyclic. 

We  define  this  graph  as  follows.  Let  5  be  a  sparsity  pattern  for  m  by  n  matrices;  in 
other  words,  S  is  a  list  of  the  entries  permitted  to  be  nonzero.  Let  G{S)  be  a  bipartite  graph 
with  one  group  of  nodes  {ri, . . . , r^}  representing  the  m  rows  and  one  group  {ci, . . . , Cn} 
representing  the  n  columns.  There  is  an  edge  from  r,-  to  Cj  if  and  only  if  Aij  is  permitted 
to  be  nonzero.  (We  will  sometimes  write  G{A)  instead  of  G{S),  where  S  is  the  sparsity 
pattern  of  /I.) 

We  also  present  another  perturbation  property  of  acyclic  matrices  which  is  quite  strong: 
multiplying  any  single  matrix  entry  by  any  factor  0^0  cannot  change  any  singular  value 
by  more  than  a  factor  of  /3  (either  up  or  down). 

Sparsity  patterns  with  this  property  have  at  most  n  +  m  -  1  nonzero  entries.  There 
are  a  great  many  such  sparsity  patterns.  Let  us  consider  only  m  by  n  sparsity  patterns  S 
which  cannot  be  permuted  into  block  diagonal  form  (this  means  G{S)  is  connected).  Then 
tlie  number  of  different  such  sparsity  patterns  is  equal  to  the  number  of  spanning  trees  on 
connected  bipartite  graphs  with  m  +  n  vertices;  this  number  is  [5,  p.  38]  [3].  If 

wo  only  wish  to  count  sparsity  patterns  which  cannot  be  made  identical  by  reordering  the 
rows  and  columns,  a  very  simple  lower  bound  on  the  number  of  such  equivalence  classes  is 
In  the  square  case  n  =  m,  Stirling’s  formula  lets  us  approximate  this 
lowoi  bound  by  £^'7(27rn^),  which  grows  quickly. 

.Since  Wo  know  the  singular  values  of  these  acyclic  matrices  are  determined  to  high 
lelativo  accuracy  by  the  data,  it  makes  sense  to  try  to  compute  them  this  accurately. 
Wo  i>iesoiu  a  bisection  algorithm  which  does  this.  The  same  algorithm  can  compute  the 
eigenvalue.')  of  arbitrary  “symmetric  acyclic”  matrices  with  tiny  componentwise  relative 
accuracy.  We  define  symmetric  acyclicity  of  a  symmetric  matrix  as  follows.  Given  a  sparsity 
pal  lorn  S  of  an  n  by  n  symmetric  matrix,  we  define  a  graph  G'{S)  by  taking  n  nodes,  and 
(  onnecting  node  i  to  node  j  ^  i  if  and  only  if  the  (i,j)  entry  is  nonzero.  The  symmetric 
xparsily  pattern  5  is  called  “symmetric  acyclic”  if  the  graph  G'{S)  is  acyclic.  (We  will 
.■ioinotinie.s  write  G'(A)  instead  of  G'{S)  where  5  is  the  sparsity  pattern  of  A.)  The  algorithm 
ovalnatos  tho  inertia  of  such  a  matrix  by  doing  symmetric  Gaussian  elimination,  with  the 
order  of  elimination  determined  by  a  postorder  traversal  of  G'{S). 

In  summary,  the  well-known  attractive  properties  of  bidiagonal  matrices  B  and  symmet¬ 
ric  iridiagonal  matrices  T,  that  the  singular  values  of  B  can  be  computed  to  high  relative 
accuracy  and  the  eigenvzJues  of  T  computed  with  tiny  componentwise  relative  backward  er¬ 
ror.  have  i)een  extended  to  “acyclic”  matrices.  In  the  case  of  computing  singular  values,  we 
have  shown  that  this  extension  is  complete:  no  other  sparsity  patterns  have  this  property. 
W'e  .strongly  suspect  that  the  set  of  symmetric  acyclic  matrices  is  also  the  complete  set  of 
symmetric  matrices  whose  eigenvalues  can  be  computed  with  tiny  componentwise  relative 
backward  error  independent  of  the  values  of  the  matrix  entries. 

Other  algorithms  for  the  special  case  of  “arrow”  matrices  are  discussed  in  [1,2,15,22]. 
I'his  work  generalizes  the  adaptations  of  bisection  to  arrow  matrices,  and  is  almost  certainly 
more  stable  than  the  QR  based  schemes. 

The  rest  of  this  paper  is  organized  as  follows.  Section  2  states  the  perturbation  theorem 
for  the  singular  values  of  acyclic  matrices,  and  section  3  proves  it.  Section  4  shows  how 
to  compute  eigenvalues  of  symmetric  acyclic  matrices  with  tiny  componentwise  relative 
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Ijackward  error,  and  applies  this  to  compute  the  singular  values  of  acyclic  matrices  to  high 
lelative  accuracy.  Section  5  give  some  examples  of  matrices  with  acyclic  sparsity  patterns. 
Section  6  discusses  algorithms  and  open  problems. 

2  Statement  of  Perturbation  Theorem  for  Singular  Values 

Ill  this  section  we  define  three  properties  of  sparsity  patterns  of  matrices,  one  about  graph 
tlieory  and  two  about  perturbation  theory.  Our  main  result,  which  we  prove  in  the  next 
section,  is  that  these  properties  are  equivalent. 

Let  A  be  an  m  by  n  matrix  with  a  fixed  sparsity  pattern  5. 

Property  1  .  6'(5)  is  acyclic. 

Property  2  .  Given  sparsity  pattern  S,  there  exist  positive  constants  to  and  with  the 
following  property.  Let  A  be  any  matrix  with  sparsity  pattern  S,  and  Aij  any  nonzero  entry. 
C'lioouf  any  |f|  <  to,  and  let  A'  =  A  except  for  =  A,j(H-£).  Then  for  all  singular  values 

(1  -  Clf|W(^)  <  Ok{A')  <  (1  +  Cl£lK(A) 

III  other  words  a  sufficiently  small  relative  perturbation  e  in  any  single  matrix  entry  cannot 
caiisi  a  relative  perturbation  greater  than  (e  in  any  singular  value. 

If  p  entries  of  A  are  simultaneously  perturbed,  Property  2  can  be  applied  p  times  to  show 
no  singular  value  can  change  by  afactor  outside  the  interval  from  (1-^|€|)p  =  l-pCkl“0(€^) 
to  ( 1  +  (,‘|f  1  )*'  =  1  +  pCkI  -  O(e^).  Property  2  seems  rather  weak,  since  it  imposes  no  bounds 
on  cy  noi  C  Still,  since  cq  and  (  are  independent  of  the  matrix  entries,  it  is  actually 
demanding  quite  a  bit  of  S.  The  last  property  is  even  stronger: 

Property  3  .  Given  sparsity  pattern  S,  let  A  be  any  matrix  with  this  sparsity,  and  Aij  any 
nonzero  entry.  Let  /?  be  any  nonzero  constant.  Let  A'  —  A  except  for  A'^j  =  PAij.  Then  for 
all  ■'<ingular  values  (7k(A') 

min(|/3|,|/?"^|)<rfe(A)  <  <Tk{A')  <  max(|/?|, |/?"^|)cr*(A) 

Property  3  is  much  stronger  than  Property  2  because  it  imposes  no  limit  to  on  the  size 
of  the  relative  perturbation,  and  because  it  asserts  C  =  1*  i-®-  the  relative  change  in 
tile  singular  values  cannot  exceed  the  relative  change  in  the  single  perturbed  matrix  entry. 
In  the  case  of  simultaneous  small  relative  perturbations  of  size  at  most  /?  =  l  +  «  in  p  entries 
of  A.  Property  3  implies  that  no  singular  value  can  change  by  a  factor  outside  the  interval 
from  (1  -  lf|)P  =  1  -  p|£|  +  to  (1  -  Icj)”'’  =  1  +  pjcj  +  (c^).  Since  the  maximum  number 
of  nonzeros  is  m  +  n  -  1,  this  relative  perturbation  is  bounded  by  (m  +  n  —  l)le|  +  O(e^). 
Our  main  result  is 

Theorem  1  Properties  1,  2  and  3  of  a  sparsity  pattern  S  are  equivalent. 
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Figure  1:  Computing  Dr  and  Dc 


if  q  is  a  row  node  then 
suppose  q  =  r, 
if  r,  is  the  root  then 
Dr,ii  =  1 

else 

suppose  Cj  is  the  parent  of  q 
Dt,u  =  ^/i-^ijDcJj) 
end 

else  {q  must  be  a  column  node)  then 
suppose  q  =  Cj 
if  Cj  is  the  root  then 

^c,jj  =  1 

else 

suppose  r,  is  the  parent  of  q 

^C,JJ  ~ 

end 
end  if 


3  Proof  of  Perturbation  Theorem  for  Singular  Values 

I'lie  proof  of  equivalence  will  consist  of  the  following  steps.  We  already  know  that  Property  3 
implies  Property  2,  so  it  will  suffice  to  prove  Property  1  implies  Property  3,  and  Property  2 
implies  Property  1. 

Lemma  1  Let  A  have  sparsity  pattern  S,  and  suppose  G{S)  is  acyclic.  Then  there  are  di- 
ayoiial  Diat  rices  Dr  and  Dc  such  that  each  entry  of  Dr  ADc  is  either  0  or  1.  Each  diagonal 
(lUry  of  Dr  or  Dc  is  a  quotient  of  monomials  in  the  entries  of  A.  In  each  monomial  each 
distinct  factor  Aij  which  appears  has  unit  exponent.  Each  Aij  can  appear  only  in  in  numer¬ 
ators  of  entries  of  Dr  and  denominators  of  entries  of  Dc,  or  vice  versa,  in  denominators 
of  entries  of  Dr  and  numerators  of  entries  of  Dc- 

Proof  Since  G(S)  is  acyclic,  it  is  a  forest  of  trees.  We  may  consider  each  tree  indepen¬ 
dently.  We  traverse  each  tree  via  depth  first  search,  and  execute  the  program  in  Figure  1 
when  first  visiting  node  q. 

The  depth  first  search  visits  each  node  once.  Since  the  graph  is  bipartite,  row  nodes  and 
column  nodes  alternate,  so  the  parent  of  a  row  node  is  a  column  node  and  vice  versa.  Since 
cacli  node  is  visited  once,  the  above  program  is  executed  once  for  each  edge  in  the  tree,  i.e. 
once  for  each  nonzero  entry  Aij,  corresponding  to  the  edge  connecting  nodes  r,-  and  Cj.  Thus 
each  Dr.u  and  Dc,jj  is  set  exactly  once.  Since  the  i,j  entry  of  Dr  ADc  is  DT,iiAijDcj},  we 
see  immediately  from  the  way  Dr^u  amd  Dc,jj  are  defined  that  this  quantity  is  1  if  Aij  ^  0 
(and  0  otherwise).  Since  each  Aij  is  used  once  during  the  graph  traversal,  each  Dr, a  and 
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must  be  be  a  quotient  of  monomiaJs.  If  Aij  is  first  used  in  Dr,ii,  then  the  formulas  in 
the  above  program  and  the  fact  the  row  and  column  nodes  alternate  mean  that  Aij  will  only 
appeal  in  denominators  of  entries  of  Dr  and  numerators  of  entries  of  Dg.  Alternatively,  if 
.-l,j  is  first  used  in  then  Aij  will  only  appear  in  denominators  of  entries  of  Dg  and 

numerators  of  entries  of  Dr.  0 

The  rest  of  the  proof  mimics  that  of  [4,  Thm.  1].  Let  E  be  the  matrix  of  ones  and  zeros 
with  sparsity  S,  so  that  DrADg  =  E.  Write  Dr  =  ^rlDrl  where  \Dr\  is  the  matrix  of  absolute 
values  of  D, .  and  Sr  is  a  diagonal  matrix  with  |5r|  =  I-  Similarly  write  Dg  =  5clZ)c|.  Then 

A  =  d;^ed:^  =  =  s;^\A\s:^ 


so  that  .4  is  related  to  |/1|  be  pre-  and  postmultiplication  by  diagonal  orthogonal  matrices. 
In  particular,  .4  and  \A\  have  the  same  singular  values.  We  wiD  henceforth  assume  without 
loss  of  generality  that  A  is  nonnegative  and  so  Dr  and  Dg  are  also  nonnegative. 

It  is  known  that  the  singular  values  of  A  are  the  same  as  the  positive  eigenvalues  of  the 
pencil 


0  A 
0 


-XI 


which  are  in  turn  the  same  as  the  positive  eigenvalues  of  the  equivalent  symmetric  definite 
pencil 


I),  0 

« 


0  A 
0 


Dr  0 
0  Dg 


0  E 
E'^  0 


-A- 


P?  0 

0  Dl 


F-XD^ 


.Now  suppose  we  perturb  A  by  changing  nonzero  entry  Aij  to  resulting  in  the 
perturbed  matrix  A' .  Apply  the  algorithm  in  Lemma  1  to  compute  a  new  P'  and  P'.  Since 
by  Lemma  1  the  entries  of  P'  and  P^  are  quotients  of  monomials  where  each  independent 
lactor  appears  at  most  once,  each  entry  P' must  equal  either  Dr,kk,  PDr.kk  or  l3~^Dr,kk- 
An  analogous  statement  about  P' and  Pc.jtfc  is  true.  Since  a  factor  Aij  must  appear  either 
in  numerators  of  Dr  and  denominators  of  Dg,  or  in  denominators  of  Dg  and  numerators  of 
Dr  -  we  have  two  cases: 


1.  Either  D'^j^^.  =  Dr,kk  or  P^  ^^fc  =  0Dr,kk,  and  either  P^  ^^^  =  ^c,kk  or  P' =  /?”^Pc,jtfc. 

2.  Either  P;  =  Dr,kk  or  P;  =  p-'^Dr,kk,  and  either  P^  =  Dg,kk  or  P^^  =  /?Pc,aa- 

Note  we  may  multiply  Dr  by  any  nonzero  7  and  divide  Dg  by  7  without  changing  the 
fact  that  DrADg  =  E.  Corresponding  to  the  above  two  cases,  we 

1.  divide  Dr  by  and  multiply  Dg  by  1/31^^^,  or 

2.  divide  Dg  by  and  multiply  Dr  by 

The  end  result  will  be  P'  and  P'  matrices  each  of  whose  entries  differs  from  the  corre¬ 
sponding  entry  of  Dr  and  Dg  by  factors  of  In  particular,  this  implies 
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loi  any  nonzero  vector  i.  Let  D'  =  diag  (DjjDj)  as  we  above  defined  D  =  diag  (X>i,Z?2). 


I’Jieii 


wr  < 


y^D'^y 

y'^D'^y 


<m 


for  any  nonzero  vector  y.  We  may  now  apply  [4,  Lemma  2]  to  conclude  that 


=  min 
8*= 


max 
a;  e  S* 

11^112  =  1 


x^Fx 

D^x 


and 

ak{A')  =  min  max 

s''  I  e  s'' 

||X||2  =  1 

svliere  tlie  minima  are  over  all  k  +  max(n,m)  dimensional  subspaces  S^,  can  differ  by  no 
more  than  a  factor  of  (i.  This  proves  that  Property  1  implies  Property  3. 


x^Fx 
x^D'^x  ’ 


Lemma  2  Let  A  have  sparsity  pattern  S,  and  let  all  its  nonzero  entries  be  independent 
iiuleienuinaies.  Then  G{S)  is  acyclic  if  and  only  if  all  minors  of  A  are  either  0  or  mono- 

minis. 


Proof  We  begin  by  noting  that  to  each  term  in  the  determinant  of  an  s  by  s  square 
matrix  .M  corresponds  a  unique  perfect  matching  in  graph  G(M).  This  is  because  each 
(01  in  in  the  determinant  corresponds  to  a  choice  of  s  '“ntries  of  M  located  in  disjoint  rows 
and  columns,  and  each  such  choice  of  s  entries  selects  a  perfect  match  in  G{M). 

Now  suppose  a  square  submatrix  M  of  A  has  at  least  two  terms  in  its  determinant, 
riiese  correspond  to  two  different  perfect  matchings.  Take  the  symmetric  difference  of  the 
edges  in  these  matchings.  This  symmetric  difference  forms  a  cycle,  which  we  get  by  following 
edges  of  tlie  tisu  matchings  in  alternation.  Thus  G{M)  contains  a  cycle,  and  so  must  G{A) 
since  it  includes  G(M). 

Now  suppose  G(v4)  contains  a  cycle.  Assume  without  loss  of  generality  that  it  is  a  simple 
cycle,  i.e.  it  is  connected  and  visits  each  node  once.  Let  M  by  the  corresponding  square 
submatrix.  This  cycle  determines  two  perfect  matchings  in  G(M),  consisting  of  alternate 
edges  of  the  cycle.  This  means  det(M)  has  at  least  two  terms.  □ 

To  prove  that  Property  2  implies  Property  1,  we  will  show  the  contrapositive.  So  assume 
6’(.d )  contains  a  cycle,  and  let  M  be  an  s  by  s  submatrix  whose  determinant  has  at  least 
2  terms.  This  means  we  may  choose  all  the  entries  of  M  to  be  nonzero  but  such  that 
M  is  exactly  singular.  Thus  its  singular  values  include  at  least  one  which  is  exactly  zero. 
Scale  M  so  that  its  entry  of  smallest  absolute  value  is  1,  and  let  a  =  |lAf||2  >  1.  Now  let 
A{Al.  If)  denote  the  matrix  with  sparsity  5,  submatrix  M,  and  other  nonzero  entries  equal 
Id  If.  Tlieu  A(M,Q)  will  have  at  least  min(m,n)  --  s  +  1  zero  singular  values,  min(m,n)  -  s 
from  the  zero  rows  and  columns  outside  M,  and  1  from  the  singularity  of  Af .  By  standard 
perturbation  theory  A{M,Tf)  will  have  at  least  rain(m,  n)  -  s  +  1  singular  values  no  larger 
tiiaii  mini.  Now  change  a  smallest  entry  of  M  from  1  to  1  +  x  to  get  M*;  thus  x  is  also  the 
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relative  change  in  this  entry.  Then  |det(A'/i)|  >  i,  and  so  OminiMx)  >  +  This 

means  cTj(.4( A/j, t}))  >  1i|/(<t  +  -  mnr],  whereas  (Ta{A{M,T]))  <  mn-q.  Thus 

<7,(A(A/r,T?))  ^  (<7>r)»4T  -  rnnq  ^  x _ ^ 

(Ts{A{M,q))  ~  mnq  mnq^a  + 

li  Property  1  held,  then  we  would  be  able  to  find  cq  >  0  and  C  >  0  such  that  for  all 
0  <  .1  <  lo  and  q  >  0  the  following  inequality  would  hold: 


mnq(a  + 

Since  we  can  make  q  as  small  as  we  like,  this  inequality  cannot  hold  for  any  finite  C-  Thus 
Propei  ty  2  cannot  hold.  This  completes  the  proof  that  Property  2  implies  Property  1,  and 
Ml  also  completes  the  proof  of  Theorem  1. 

4  A  bisection  algorithm  for  computing  eigenvalues  with 
tiny  backward  error 

l.et  :  denote  the  machine  precision.  We  will  assume  the  usual  model  of  floating  point 
error.  //(«  ,  6)  =  (a  0  6)(1  +  i)  with  <  £a,,  and  assume  neither  underflow  nor  overflow 
occur.  (Of  course,  a  practical  algorithm  would  need  to  account  for  overflow.  This  can 
lie  done  analogously  to  the  way  overflow  is  accounted  for  in  standard  tridiagonal  bisection 

In  this  section  we  will  show  how  to  compute  the  eigenvalues  of  a  symmetric  acyclic 
maiii.x  7  with  tiny  componentwise  relative  backward  error.  Our  main  result  is 

Theorem  2  fhe  algorithm  in  Figure  2  computes  count(T,x),  the  number  of  eigenvalues  of 
I  /c.s,^  llutii  .r.  with  a  backward  error  6T  with  the  following  properties: 

\yj',,\  <  (1..5c  +  2.5)£Ml7’.ji  u-’^en  i  j. 

|AT„1  <  (2i’  + 2)5m1i1.  : 

Ht  ri  (•<;/-  1  is  the  maximum  degree  of  any  node  in  the  graph  of  T.  In  other  words,  the 
i  oiiiputtd  count(T.x)  is  the  exact  value  of  count(J’  +  6T,x)  where  ST  is  bounded  as  above. 

This  is  essentially  identical  to  the  standard  error  analysis  of  Sturm  sequence  evaluation 
lor  symmetric  tridiagonal  matrices  [9,  Sec.  6]  [13]  (this  is  stronger  than  the  result  in  [20,  p. 
.1031). 

Our  algorithm  simply  performs  symmetric  Gaussian  elimination  on  T  —  xl\  P{T  - 
xl)P^  =  LDLJ  where  P  is  a  permutation  matrix,  L  is  unit  lower  triangular  and  D  is 
diagonal.  Then  count(T', x)  is  simply  the  number  of  negative  diagonal  entries  of  D,  by 
.Sylvester  s  Inertia  Theorem  [16].  The  order  of  elimination  is  the  same  as  a  postorder 
t  ra\ersal  of  the  nodes  of  the  acyclic  graph.  Since  leaves,  which  have  degree  1,  are  eliminated 
first,  there  is  no  fill-in  during  the  elimination,  and  all  off-diagonal  entries  Lij  of  L  can  be 
computed  by  simply  dividing  Lij  —  TijJDjj. 
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Figure  2:  Computing  count(r,x) 


call  ci\l(i,d,s,x)  where  i  is  any  node  1  <  i  <  n 
leluni  count(T, i)  =  a 

procedure  cnt(i,d,  s,i) 

/*  i  and  x  are  input  parameters,  d  and  s  are  output  parameters  */ 
d  =  T„  -  X 

•s  =  0 

for  all  children  j  of  i  do 
call  cnl(j,d' ,s',x) 
d  =  d-  Tf^/d' 

,s  =  -f-  s' 

end  for 

if  d  <  0,  then  a  =  a  +  1 
return  d  and  s 
end  procedure 


Wo  a^^^une  the  graph  G'{S)  is  connected,  since  otherwise  the  matrix  can  be  reordered  to 
1)0  block  diagonal  (one  diagonal  block  per  connected  component  of  G'(5)),  and  the  inertia 
of  each  diagonal  block  can  be  computed  separately.  The  algorithm  cnt(t,d,s,i)  in  Figure  2 
assumes  the  matrix  is  stored  in  graph  form.  Subroutine  cnt(t,d,s,i)  does  a  postorder 
travel  sal  of  the  acyclic  graph  G'(5),  and  may  be  called  starting  at  any  node  1  <  t  <  n.  In 
addition  to  i,  x  is  an  input  parameter.  The  variables  d  and  s  are  output  parameters;  on 
return  .s  ir,  the  desired  value  of  count(T,i). 

I'o  pro\e  Theorem  2,  we  will  exploit  the  acyclicity  of  T  to  show  that  each  computed 
c|uaiitit.\-  and  original  entry  of  T  is  used  (directly)  just  once  during  the  entire  computation, 
and  then  use  this  to  “push"  the  rounding  error  back  to  the  original  data. 

We  see  that  each  entry  of  T  is  used  just  once  as  follows.  Tu  is  only  used  when  visiting 
node  i.  and  Tij  is  used  only  once,  when  visiting  i  if  j  is  a  child  of »  or  when  visiting  j  if  i  is 
a  child  of  j  in  the  postorder  traversal  tree. 

.Now  denote  the  d  computed  when  visiting  node  i  by  d,.  The  floating  point  operations 
performed  while  visiting  node  i  are  then 


d,  =  fl 


Ti,  -  I  - 


E 

all  children 
j  of  i 


(4.1) 


I'o  analyze  this  formula,  we  will  let  subscripted  e’s  denote  independent  quantities 
l)ounded  in  absolute  value  by  s*#-  We  will  also  make  standard  approximations  like 
il  +f,)*'(l  +r2)*‘  =  1  +  2£3. 
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Siluo  we  do  not  know  the  number  of  terms  or  the  order  of  the  sum  in  equation  (4.1),  we 
will  make  the  worst  case  assumption  that  there  are  v  <n-  \  terms  where  v  is  the  maximum 
(lef>iee  of  any  node  in  the  graph  G'{S).  This  leads  to 


fl,  —  il  +  {v  +  l)Sia)Tii  —  (1  +  (n  +  l)^,*)^  —  ^  (!  +  («  +  (4-2) 

all  chi’a.'en  ^ 

j  of  i 


71  77T  n~  ~  3-+  (2r  +  2)£,cX  -  2-  - 1 

allcMdrcn 
j  of  i 

l.ci  ;  ,,,  he  tlie  roundoff  error  corresponding  to  Sia  committed  when  computing  dj.  Then 


1  +  (  ('  +  1  )5,(i 


=  T„  -  X  +  (2n  +  2)£icX 


all  children 
j  of  i 


((l  +  (1.5t>  +  2.5)£.j»)r.;)^ 
dj/{l-\-{v+l)£ja) 


\\  finalh. 


tl',  =  -  X  +  {2t'  +  2)e,^i  -  ^ 

all  children 
j  of  i 


((l  +  (1.5v  +  2.5)£.j»)r,j)^ 


where  d'  =  (/,/( 1  +  fia).  Equation  (4.5)  tells  us  that  the  d[  are  the  exact  diagonal  entries  of 
I)  ill  P('r  +  6T  -  xI)P^  =  LDL^ .  Since  they  obviously  have  the  same  signs  as  the  d,,  this 
|)rove^  'riieorem  2. 

Tlie  proof  depends  strongly  on  there  not  being  any  fill-in  and  on  each  off  diagonal  entry 
being  computable  by  a  single  division.  Since  these  properties  hold  if  and  only  if  the  graph 
(i'{  7  )  is  symmetric  acyclic,  we  strongly  suspect  that  this  is  the  only  class  of  matrices  whose 
eigenvalues  can  always  be  computed  with  tiny  componentwise  relative  backward  error. 

We  now  apply  Theorem  2  to  compute  singulcir  values  of  acyclic  matrices  to  high  relative 
accuracy.  So  suppose  B  is  a  matrix  whose  graph  G{B)  is  acyclic.  Consider  the  symmetric 


matrix 

0  B 
B^  0 


It  is  wei  known  that  the  positive  eigenvalues  of  A  are  the  singular  values  of  B.  It  is 
also  immediate  that  the  graph  G'{A)  =  G{B).  Therefore  B  is  acyclic  if  and  only  if  A  is 
s\  inmetric  acyclic,  so  we  can  apply  the  above  algorithm  to  compute  all  B'&  singular  values 
to  liigh  relative  accuracy. 

One  other  algorithm  is  worth  mentioning.  If  4  is  symmetric  positive  definite  and  sym¬ 
metric  acyclic,  then  its  Cholesky  factor  L  is  acyclic,  has  the  “lower  half”  of  the  sparsity 
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pattern  of  and  may  be  computed  by  using  algorithm  cut.  It  may  occasionally  be  more 
accurate  to  compute  yl’s  eigenvalues  by  first  computing  L,  computing  its  singular  values  by 
bisection,  and  then  squaring  the  singular  values  to  get  j4’6  eigenvalues  [4].  This  is  the  case, 
for  example,  for  the  tridiagonal  matrix  with  2’s  on  the  diagonal  and  I’s  on  the  ofT-diagonal. 

5  Examples 

\V(>  give  \arious  examples  of  acyclic  sparsity  patterns,  beginning  with  acyclic  G(S).  Given 
any  acyclic  sparsity  pattern,  others  can  be  generated  either  by  permuting  rows  and/or 
cohimn.s,  or  by  adding  more  zeros.  Since  all  square  acyclic  matrices  have  monomial  (or 
zero)  determinants,  this  means  we  can  permute  them  to  be  upper  triangular.  In  addition 
to  bidiagona]  matrices,  some  other  examples  are 


_ 

'  X  X 

X 

X 

X  X 

X  X 

X  X 

X  X 

X 

and 

X  X 

X  X 

X  X 

X  . 

.1  = 


I'o  get  symmetric  acyclic  matrices  A,  one  can  always  take  an  acyclic  B  and  set 
^  I  -  A/.  Some  other  examples  are 


0  J 


6  Algorithms  and  Open  Problems 

In  [8]  a  perturbation  theorem  for  singular  vectors  of  bidiagonal  matrices  is  proven,  which 
•shows  that  the  appropriate  condition  number  for  the  t-th  singular  vector  is  the  reciprocal 
of  the  relative  difference  between  the  i-th  singular  value  and  next  closest  one.  It  wotild  be 
interesting  to  extend  this  to  the  acyclic  case. 

Given  the  perturbation  theory,  it  would  be  nice  to  compute  the  singular  vectors  as 
accurately  as  they  deserve.  A  natural  candidate  is  inverse  iteration,  but  even  in  the  simple 
case  of  symmetric  tridiagonal  matrices,  open  problems  remain.  In  particular  there  is  no 
absolute  guarantee  that  the  computed  eigenvectors  are  orthogonal,  although  in  practice  the 
algorithm  can  be  made  quite  robust  [11]. 

In  the  ■’extreme”  cases  of  tridiagonal  and  arrow  matrices,  we  know  how  to  compute  the 
inertia  in  Oflogn)  time,  using  the  so-called  parallel-prefix  algorithm  in  the  tridiagonal  case 
[17.19] .  and  more  simply  in  the  arrow  case.  The  stability  in  the  tridiagonal  case  is  unknown. 
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but  ill  practice  it  appears  to  be  stable.  We  can  extend  this  to  the  general  symmetric  acylic 
case  in  two  ways.  First,  the  tree  describing  the  expression  whose  final  value  is  di  htis  at  most 
»  leaves.  From  [6]  we  know  any  such  expression  tree  can  be  evaluated  in  at  most  4log2n 
parallel  steps,  although  stability  may  be  lost.  Another  approach,  which  includes  parallel 
prefix  and  the  algorithm  in  [15]  as  special  cases,  is  based  on  [14].  The  idea  is  to  simply 
evzduate  the  tree  greedily,  summing  k  leaves  of  a  single  node  in  O(log2  k)  steps  whenever 
possible,  and  collapsing  a  chain  of  k  nodes  into  a  single  node  via  parallel  prefix  in  0(log2  k) 
steps  whenever  possible.  If  we  could  understand  the  numerical  stability  of  parallel  prefix, 
we  could  probably  analyze  this  more  general  scheme  as  well. 

Divide  and  conquer  [7,10,18,12]  has  been  widely  used  for  the  tridiagonal  eigenproblem 
and  bidiagonaJ  singular  value  decomposition.  This  can  be  straightforwardly  extended  to 
the  acyclic  case.  In  terms  of  the  tree,  just  remove  the  root  by  a  “rank  one  tearing”,  solve 
the  independent  child  subtrees  recursively  and  in  parallel,  and  merge  the  results  by  solving 
the  secular  equation  [21].  Any  node  can  be  the  root,  and  to  be  efficient  it  is  important  that 
no  subtree  be  large.  In  the  tridiagonal  case,  there  are  always  two  subtrees  of  nearly  equal 
size.  In  a  general  tree  one  can  only  make  sure  that  no  subtree  has  more  than  half  the  nodes 
of  the  original  tree  (this  is  easily  done  in  0(n)  time  via  depth  first  search). 

QR  does  not  appear  to  extend  beyond  the  tridiagonal  case.  The  case  of  arrow  matrices 
was  analyzed  in  [2],  where  it  was  shown  that  no  QR  algorithm  could  exist.  A  simpler  proof 
arises  from  noting  that  two  steps  of  LL^  is  equivalent  to  one  step  of  QR  in  the  positive 
definite  case,  and  so  the  question  is  whether  the  sparsity  pattern  of  To  =  is  the  same 
as  that  of  Ti  =  L\  this  is  easily  seen  to  include  only  tridiagonal  To. 
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